!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
!
! An interface of bitmap image utility for real4 and 
! real8 data.
! Users be expected to 'use buty' instead of 'use butyD' 
! or the other butys. 
!
! TTD:
! 1) catch allocation failure
! 2) error code and error handling routine
!
! note:
! 1) some routines use one dimenstional accessing
!    to 3D array, which may cause some problems
!    when compiling with boundary check option
!    turned on.
!
module buty
  use butyD
  use butyR4
  use butyR8
  use butyI4
  use butyI2
  use butyI1
  implicit none
  interface mask_img
     module procedure mask_img_r4
     module procedure mask_img_r8
     module procedure mask_img_i4
     module procedure mask_img_i2
     module procedure mask_img_i1
  end interface
  interface median_filter
     module procedure median_filter_r8
     module procedure median_filter_r4
     module procedure median_filter_i4
     module procedure median_filter_i2
     module procedure median_filter_i1
  end interface  
  interface get_minmax_img
     subroutine get_minmax_img_i1(n,img,pmin,pmax)
       integer*4,intent(in)::n
       integer*1,intent(in)::img
       real*8,intent(out)::pmin,pmax
     end subroutine get_minmax_img_i1
     subroutine get_minmax_img_i2(n,img,pmin,pmax)
       integer*4,intent(in)::n
       integer*2,intent(in)::img
       real*8,intent(out)::pmin,pmax
     end subroutine get_minmax_img_i2
     subroutine get_minmax_img_i4(n,img,pmin,pmax)
       integer*4,intent(in)::n
       integer*4,intent(in)::img
       real*8,intent(out)::pmin,pmax
     end subroutine get_minmax_img_i4
     subroutine get_minmax_img_r4(n,img,pmin,pmax)
       integer*4,intent(in)::n
       real*4,intent(in)::img
       real*8,intent(out)::pmin,pmax
     end subroutine get_minmax_img_r4
     subroutine get_minmax_img_r8(n,img,pmin,pmax)
       integer*4,intent(in)::n
       real*8,intent(in)::img
       real*8,intent(out)::pmin,pmax
     end subroutine get_minmax_img_r8
  end interface
  interface get_avevar_img
     module procedure get_avevar_img_r4
     module procedure get_avevar_img_r8
     module procedure get_avevar_img_i4
     module procedure get_avevar_img_i2
     module procedure get_avevar_img_i1
  end interface
  interface get_hist
     module procedure get_hist_r4
     module procedure get_hist_r8
     module procedure get_hist_i4
     module procedure get_hist_i2
     module procedure get_hist_i1
  end interface
  interface kmclusterize_img
     module procedure kmclusterize_img_r4
     module procedure kmclusterize_img_r8
     module procedure kmclusterize_img_i4
     module procedure kmclusterize_img_i2
     module procedure kmclusterize_img_i1
  end interface
  interface get_2dhist
     module procedure get_2dhist_r4
     module procedure get_2dhist_r8
     module procedure get_2dhist_i4
     module procedure get_2dhist_i2
     module procedure get_2dhist_i1
  end interface
  interface remaphist_img
     module procedure remaphist_img_r4
     module procedure remaphist_img_r8
     module procedure remaphist_img_i4
     module procedure remaphist_img_i2
     module procedure remaphist_img_i1
  end interface
  interface eqhist_img
     module procedure eqhist_img_r4
     module procedure eqhist_img_r8
     module procedure eqhist_img_i4
     module procedure eqhist_img_i2
     module procedure eqhist_img_i1
  end interface
  interface contrast_norm_img
     module procedure contrast_norm_img_r4
     module procedure contrast_norm_img_r8
     module procedure contrast_norm_img_i4
     module procedure contrast_norm_img_i2
     module procedure contrast_norm_img_i1
  end interface
  interface infopot_img
     module procedure infopot_img_r4
     module procedure infopot_img_r8
     module procedure infopot_img_i4
     module procedure infopot_img_i2
     module procedure infopot_img_i1
  end interface
  interface adapbin_filter
     module procedure adapbin_filter_r4
     module procedure adapbin_filter_r8
     module procedure adapbin_filter_i4
     module procedure adapbin_filter_i2
     module procedure adapbin_filter_i1
  end interface
  interface filter_img
     module procedure filter_img_i1
     module procedure filter_img_i2
     module procedure filter_img_i4
     module procedure filter_img_r4
     module procedure filter_img_r8
  end interface
  interface get_H
     module procedure get_H_i1
     module procedure get_H_i2
     module procedure get_H_i4
     module procedure get_H_r4
     module procedure get_H_r8
  end interface
  integer*4,parameter:: SIDF_SET_MINMAX = Z'00000001'
  integer*4,parameter:: SIDF_SET_AVEVAR = Z'00000002'
  integer*4,parameter:: SIDF_SET_NSAB   = Z'00000004'
  integer*4,parameter:: SIDF_ADJ_LOUP   = Z'80000000'
  integer*4,parameter:: SIDF_SET_ALL    = Z'0FFFFFFF'
contains
  !------------------------------------------------------!
  subroutine set_image_data(x,pvar_ns_ratio,flag)
    type(image_data),intent(inout)::x
    real*8,intent(in),optional::pvar_ns_ratio
    integer*4,intent(in),optional::flag
    integer*4 f
    if(present(flag)) then
       f=flag
    else
       f=SIDF_SET_ALL
    endif
    x%xmin = -x%d*dble(x%n)/2
    x%xmax =  x%d*dble(x%n)/2
    x%offset=dble(x%n+1)*x%d/2
    if(iand(f,SIDF_ADJ_LOUP)/=0) then
       x%pup=huge(1.d0)
       x%plo=-huge(1.d0)
    end if
    if(iand(f,SIDF_SET_MINMAX)/=0) then
       select case(x%datatype)
       case(DT_IKIND_1)
          call get_minmax_img(product(x%n),x%img_i1(1,1,1),x%pmin,x%pmax)
       case(DT_IKIND_2)
          call get_minmax_img(product(x%n),x%img_i2(1,1,1),x%pmin,x%pmax)
       case(DT_IKIND_4)
          call get_minmax_img(product(x%n),x%img_i4(1,1,1),x%pmin,x%pmax)
       case(DT_RKIND_4)
          call get_minmax_img(product(x%n),x%img_r4(1,1,1),x%pmin,x%pmax)
       case(DT_RKIND_8)
          call get_minmax_img(product(x%n),x%img_r8(1,1,1),x%pmin,x%pmax)
       end select
       x%pup=min(x%pup,x%pmax)
       x%plo=max(x%plo,x%pmin)
    end if
    if(iand(f,SIDF_SET_AVEVAR)/=0) then
       select case(x%datatype)
       case(DT_IKIND_1)
          call get_avevar_img(x%n,x%img_i1,x%plo,x%pup,x%pave,x%pvar)
       case(DT_IKIND_2)
          call get_avevar_img(x%n,x%img_i2,x%plo,x%pup,x%pave,x%pvar)
       case(DT_IKIND_4)
          call get_avevar_img(x%n,x%img_i4,x%plo,x%pup,x%pave,x%pvar)
       case(DT_RKIND_4)
          call get_avevar_img(x%n,x%img_r4,x%plo,x%pup,x%pave,x%pvar)
       case(DT_RKIND_8)
          call get_avevar_img(x%n,x%img_r8,x%plo,x%pup,x%pave,x%pvar)
       end select
    endif
    if(iand(f,SIDF_ADJ_LOUP)/=0) then
       if(x%pmin>PVAL_ALLOWABLE_NEGATIVE) then
          x%plo=0
       else
          x%plo=max(x%pmin,-1000.d0) ! CT
       end if
       x%plo=max(x%pave-x%pvar*1.5,x%plo)
       x%pup=min(x%pave+x%pvar*1.5,x%pup)
!!$       if(x%pmin<PVAL_ALLOWABLE_NEGATIVE) then
!!$          x%plo=-1000.
!!$          x%pup= 1000.
!!$       else
!!$          x%plo=max(0.,x%pmin)
!!$          x%pup=x%pmax
!!$          adj_up=.true.
!!$       endif
    end if
!!$    if(adj_up) x%pup=min(x%pave+x%pvar,x%pmax) ! higher value is less significant
    if(iand(f,SIDF_SET_NSAB)/=0) then
!!$       if(present(pvar_ns_ratio)) then
!!$          x%ns=x%pvar*pvar_ns_ratio
!!$       else
!!$          x%ns=nint(x%pup-x%plo)+1
!!$       endif
       x%ns=x%pup-x%plo+1
       if(x%cf/=0) x%ns=x%ns/x%cf
       x%ns=max(GL_MIN,min(GL_MAX,x%ns))
       x%a=dble(x%ns-1)/(x%pup-x%plo)
       x%b=(x%plo*dble(x%ns)-x%pup)/(x%plo-x%pup)
    endif
  end subroutine set_image_data
  !======================================================
  subroutine mv_img(rtp,inttype,in,out)
    type(rigid_transform_param),intent(in)::rtp
    integer*4,intent(in)::inttype
    type(image_data),intent(in)::in
    type(image_data),intent(out)::out
    select case(in%datatype)
    case(DT_IKIND_1)
       call mv_img_i1(rtp,inttype,in,out)
    case(DT_IKIND_2)
       call mv_img_i2(rtp,inttype,in,out)
    case(DT_IKIND_4)
       call mv_img_i4(rtp,inttype,in,out)
    case(DT_RKIND_4)
       call mv_img_r4(rtp,inttype,in,out)
    case(DT_RKIND_8)
       call mv_img_r8(rtp,inttype,in,out)
    end select
  end subroutine mv_img
  !------------------------------------------------------!
  subroutine rebin_img(in,out)
    type(image_data),intent(in)::in
    type(image_data),intent(out)::out
    select case(in%datatype)
    case(DT_IKIND_1)
       call rebin_img_i1(in,out)
    case(DT_IKIND_2)
       call rebin_img_i2(in,out)
    case(DT_IKIND_4)
       call rebin_img_i4(in,out)
    case(DT_RKIND_4)
       call rebin_img_r4(in,out)
    case(DT_RKIND_8)
       call rebin_img_r8(in,out)
    end select
  end subroutine rebin_img
  !------------------------------------------------------!
  subroutine label2d_img(in,out,kmax)
    type(image_data),intent(in)::in
    type(image_data),intent(out)::out
    integer*4,intent(out)::kmax
    select case(in%datatype)
    case(DT_IKIND_1)
       call label2d_img_i1(in,out,kmax)
    case(DT_IKIND_2)
       call label2d_img_i2(in,out,kmax)
    case(DT_IKIND_4)
       call label2d_img_i4(in,out,kmax)
    case(DT_RKIND_4)
       call label2d_img_r4(in,out,kmax)
    case(DT_RKIND_8)
       call label2d_img_r8(in,out,kmax)
    end select
  end subroutine label2d_img
  !------------------------------------------------------!
  real*8 function get_mi_w_rtp(rtp,ref,tar)
    type(rigid_transform_param),intent(in)::rtp
    type(image_data),intent(in)::ref
    type(image_data),intent(in)::tar
    real*8 a(3,4)
    call set_affine_rigid(rtp,a)
    select case(ref%datatype)
    case(DT_IKIND_1)
       get_mi_w_rtp=get_mi_w_affine_i1(a,ref,tar)
    case(DT_IKIND_2)
       get_mi_w_rtp=get_mi_w_affine_i2(a,ref,tar)
    case(DT_IKIND_4)
       get_mi_w_rtp=get_mi_w_affine_i4(a,ref,tar)
    case(DT_RKIND_4)
       get_mi_w_rtp=get_mi_w_affine_r4(a,ref,tar)
    case(DT_RKIND_8)
       get_mi_w_rtp=get_mi_w_affine_r8(a,ref,tar)
    end select
  end function get_mi_w_rtp
  !------------------------------------------------------!
  real*8 function get_chi2_w_rtp(rtp,ref,tar)
    type(rigid_transform_param),intent(in)::rtp
    type(image_data),intent(in)::ref
    type(image_data),intent(in)::tar
    real*8 a(3,4)
    call set_affine_rigid(rtp,a)
    select case(ref%datatype)
    case(DT_IKIND_1)
       get_chi2_w_rtp=get_chi2_w_affine_i1(a,ref,tar)
    case(DT_IKIND_2)
       get_chi2_w_rtp=get_chi2_w_affine_i2(a,ref,tar)
    case(DT_IKIND_4)
       get_chi2_w_rtp=get_chi2_w_affine_i4(a,ref,tar)
    case(DT_RKIND_4)
       get_chi2_w_rtp=get_chi2_w_affine_r4(a,ref,tar)
    case(DT_RKIND_8)
       get_chi2_w_rtp=get_chi2_w_affine_r8(a,ref,tar)
    end select
  end function get_chi2_w_rtp
  !------------------------------------------------------!
  real*8 function get_mi_w_nrtp(rtp,ref,tar)
    type(rigid_transform_param),intent(in)::rtp
    type(image_data),intent(in)::ref
    type(image_data),intent(in)::tar
    real*8 a(3,4)
    call set_affine(rtp,a)
    select case(ref%datatype)
    case(DT_IKIND_1)
       get_mi_w_nrtp=get_mi_w_affine_i1(a,ref,tar)
    case(DT_IKIND_2)
       get_mi_w_nrtp=get_mi_w_affine_i2(a,ref,tar)
    case(DT_IKIND_4)
       get_mi_w_nrtp=get_mi_w_affine_i4(a,ref,tar)
    case(DT_RKIND_4)
       get_mi_w_nrtp=get_mi_w_affine_r4(a,ref,tar)
    case(DT_RKIND_8)
       get_mi_w_nrtp=get_mi_w_affine_r8(a,ref,tar)
    end select
  end function get_mi_w_nrtp
end module buty
